%% correct_ica_pvals.m
%
% [corrected_pvals scaling]=correct_ica_pvals(pvals,threshold, strength)
% Input arguments:  pvals: column vector of pvalues
%                   p-value threshold: default 0.05
%                   strength: a string reading either 'weak' (default) or
%                   'strong'. Weak correction employs a multi-step-up test.
%                   Strong employs Bonferroni correction.
%                
% Output arguments: corrected_pvals: multi-step-up corrected p-values
%                   (corrected_pvals=pvals*scaling)
%                   scaling: scaling factor for multi-step-up correction.
% HL 22.05.12

function [corrected_pvals scaling]=correct_ica_pvals(pvals,threshold,strength)

if nargin<3;
    strength='weak';
end
if nargin<2;
    warning('Threshold not specified. Set to 0.05 by default')
    threshold=0.05;
end
if nargin<1
    error('No p-values provided')
end

switch lower(strength)
    case {'weak'}
        fprintf('\nApplying weak multi-step-up correction to p-values\n ')
        [pvals_sorted order]=sort(pvals,1,'descend');
        i=1;
        while pvals_sorted(i)*length(pvals)/i>threshold
            i=i+1;
            
            if i>length(pvals_sorted)
                error('P-values too big. No FDR correction possible')
            end
            
        end
        scaling=length(pvals)/i;
        corrected_pvals=scaling*pvals;
    case {'strong'}
        fprintf('\nApplying Bonferroni correction to p-values\n ')
        corrected_pvals=pvals*numel(pvals);
        scaling=ones(numel(pvals),1)*numel(pvals);
end
end